Properties making a chaotic system a good Pseudo Random Number Generator 
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We discuss the properties making a deterministic algorithm suitable to generate a pseudo random 
sequence of numbers: high value of Kolmogorov- Sinai entropy, high-dimensionality of the parent 
dynamical system, and very large period of the generated sequence. We propose the multi dimen- 
sional Anosov symplectic (cat) map as a Pseudo Random Number Generator. We show what chaotic 
features of this map are useful for generating Pseudo Random Numbers and investigate numerically 
which of them survive in the discrete state version of the map. Testing and comparisons with other 
generators are performed. 

PACS numbers: 05.45.Jn, 05.45.Pq, 05.40.-a 



I. INTRODUCTION 

In most scientific uses of numerical computations, 
e.g. Montecarlo simulations and molecular dynamics, 
it is necessary to have a series of independent, iden- 
tically distributed (i.i.d.) continuous random variables 
x(l),x(2), . . . , x(n) with assigned single variable proba- 
bility density function (pdf) P(x(i)). Of course, it is 
enough to have i.i.d. random variables {x(t)} uniformly 
distributed in the interval [0,1], since a suitable change 
of variable y = g{x) may generate numbers {y{i)} with 
any pdf P(y). 

Let us call perfect random number generator (RNG) a 
process producing i.i.d. variables uniformly distributed 
in [0,1]. One can produce perfect RNG only using non- 
deterministic physical phenomena, e.g. the decay of ra- 
dioactive nuclei or the arrival on a detector of cosmic 
rays. 

A more practical way is to use a computer that 
produces a "random-looking" sequence of numbers, by 
means of a recursive rule. Let us call Pseudo Random 
Number Generator (PRNG) an algorithm designed to 
mimic a random sequence on a computer. This issue is 
far from being trivial; in Von Neumann words: "Anyone 
who considers arithmetical methods of producing random 
digits is, of course, in a state of sin" pj. The two un- 
avoidable problems are the following: 

a) Numerical algorithms are deterministic. 

b) They deal with discrete numbers. 

The limitations arising from these properties can be ana- 
lyzed using the language and the tools of dynamical sys- 
tems theory. In the following, we anticipate how these 
remarks translate in this framework and the main issues 



of the entropic characterization of PRNG's (these issues 
are discussed in detail in Sec. HT|) : 

a\) since the algorithm is deterministic, the Kol- 
mogorov - Sinai (KS) entropy (has) is finite. The 
sequence {x(i)} can not be "really random", i.e. 
with an infinite KS entropy, because the determin- 
istic dynamical rule constrains the outputs that are 
near in time and supplies us with a maximum of 
\og 2 (e hKS ) random bits per unit time. This limita- 
tion would be present also in a hypothetical com- 
puter able to work with real numbers. 

b\) Since any deterministic system with a finite number 
of states is periodic, any sequence produced by an 
algorithm working with discrete numbers must be 
periodic, possibly after a transient: therefore, not 
only hxs < °°, but hxs = 0. The computer- 
implemented system can be only pseudo-chaotic. 

Firstly, we consider point (a). After the seminal work 
of Lorenz and Henon 0, Q (to mention just two of the 
founders of the modern theory of chaos), it is well estab- 
lished that also deterministic systems may have a time 
evolution that appears rather "irregular" with the typi- 
cal features of genuine random processes. This evidence 
opened a debate on the possibility to distinguish between 
noisy and chaotic deterministic dynamics. Following the 
work of Takens |j], the so-called embedding techniques 
have been developed to extract qualitative and quantita- 
tive information from a time series. An initial enthusiasm 
was due to that the use of the embedding method (via 
delayed coordinates) allows, at least in principle, the de- 
termination of quantities like dimensions, hxs & n d Lya- 
punov exponents. People believed that, after determin- 
ing the KS entropy of a data sequence, one would know 
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the true nature (deterministic or stochastic) of the law 
generating the series. It is now rather clear that there 
are several limitations in the use of this technique p| : for 
instance, the number of points necessary for the phase- 
space reconstruction increases exponentially with the di- 
mension of the system gj. Thus, due to the finiteness 
of the datasets, it is not possible to perform an entropic 
analysis with an arbitrary fine resolution, i.e. to compute 
the e-entropy h(e) for very small values of e. This fact 
severely restricts the possibility to distinguish between 
signals generated by different rules, such as regular (high 
dimensional) systems, deterministic chaotic systems, and 
genuine stochastic processes. Although the above result 
may appear negative, it allows a pragmatic classifica- 
tion of the stochastic or chaotic feature of the signal, 
according to the dependence of the e-entropy on e and 
this yields some freedom in modeling systems Q. As 
a relevant example of a representation of a determinis- 
tic system in term of stochastic processes, we mention 
the fully developed turbulence @|. Turbulent systems 
are high dimensional deterministic chaotic systems and 
therefore h(e) ~ tigs f° r e % £c, where e c — > as the 
Reynolds number Re — > oo, while h(e) ~ e~ a for e > e c . 
The fact that in certain stochastic processes h(e) ~ e~ a 
can be useful for modeling purposes: for example, in the 
so-called synthetic turbulence, one introduces suitable 
multi-affine stochastic processes with the correct scaling 
properties of the fully developed turbulence. 

In this paper we want to discuss about the opposite 
strategy, i.e. to mimic noise with deterministic chaotic 
systems. Let us summarize the starting points of our ap- 
proach to use a deterministic chaotic system as a PRNG: 

1. Since in any deterministic system h(e) ~ has for 
e < e c with (In eg) ~ —h-KS, one should work with a 
very large hj^s 9] . In this way the true (determin- 
istic) nature of the PRNG becomes apparent only 
at a very high resolution. 

2. The outputs {x(i}} of a perfect RNG, when ob- 
served at resolution e, supply log 2 (l/e) oc h(e) ran- 
dom bits per iteration. In order to observe the be- 
havior h(e) ~ ln(l/e) for e > e c in a deterministic 
algorithm, it is necessary that the time correlation 
is very weak. We will discuss how this property 
may be achieved by taking as output a single vari- 
able of a high-dimensional chaotic system. 

It is not difficult to satisfy point (1), while request (2) is 
less obvious. Anosov systems 10] are natural candidates 
to fulfill it, having invariant stationary measure and very 
strong chaotic properties. 

A third point has to be added, dealing with the prob- 
lem (b) and its consequence (&i). 

Up to here we were considering the chaotic properties 
of systems with continuous phase space. Quantities like 
the hxs and the e-entropy have an asymptotic nature, i.e. 
they are related to large time behavior. However, there 
are situations where the system is, strictly speaking, non 



chaotic (hies — 0) but its features appear irregular to 
a certain extent. Such property (denoted with the term 
pseudo-chaos [Til, [H, Ellis basically due to the presence 
of long transient effects [l4| . 

As noted above, the use of a computer discretizes the 
phase space of a dynamical system, canceling (at least) 
its asymptotic chaotic properties. However, if the period 
of the realized sequence is long enough, the effects re- 
lated to points (1) and (2) reasonably survive as a chaotic 
transient. According to this observation, we add a third 
request: 

3. the period of the series generated by the computer 
(i.e. with a state-discretization of the deterministic 
system) must be very large. 

Point (3) is really tough: as far as we know, for 
a generic deterministic system with Ai discrete states, 
there are no general methods to determine a priori 
the length of the periodic orbits. A nice result, based 
on probabilistic considerations, suggests that the period 
T ~ M 1 / 2 jig, although strong fluctuations are present. 
The use of high-dimensional systems may be a natural 
solution also for this problem: calling M the number of 
states along each of the d dimensions, the typical period 
T ~ M d / 2 grows very fast by increasing M and d. 

Whatever the mechanism for producing the pseudo- 
chaotic transient, the mere fact that the sequence is pe- 
riodic implies that it is possible to obtain equidistributcd 
words only up to a length m = O(lnT). Thus, long 
time correlations among outputs of a generator can not 
be detected by the standard entropic analysis. We will 
show that a high dimensional chaotic system provides 
outputs which are not correlated even looking at time 
delays greater than m. In particular we will discuss the 
connection between correlation functions and the spec- 
tral test for random sequence 0, [l8L ITflj showing that 
the outputs of the high-dimensional cat map have zero 
ri-points correlation functions, when n is less or equal to 
the dimensions of the map. 

In Section [Hj we describe the entropic properties of 
PRNGs currently used, underlying both the mechanisms 
involved in PRNGs. In Section ITTll the algorithms used 
to test the generators are described. In Section IIVI we 
study the properties of the multi-dimensional Arnol'd's 
cat map and in Section we propose its discrete version 
as a PRNG. Section IVII is devoted to conclusions and 
perspectives. 



II. ENTROPY AND GOOD PRNG 

First of all, we briefly recall some basic notion on the e- 
entropy ^3|. Consider the variable x(t) G M d represent- 
ing the state of a <i-dimensional system, and introduce 
the new variable 

(t) = (x(t), x(t + 1), . . . ,x(t + m - 1)) G K md . (1) 
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Of course, y(™ 1 ' corresponds to a trajectory in a time in- 
terval to. Then, the phase space is partitioned in cells 
of linear size e in each of the d directions. Since the re- 
gion where a bounded trajectory evolves contains a finite 
numbers of cells, each y( m ) (t) defined in can be coded 
into a word of length m out of a finite alphabet: 

y M (t) _ wW(t) = (*(e, t), t(e, . . . , i(e, t+m-1)) 

(2) 

where i(e, t + j) labels the e-cell containing x(t + j). 
Assuming that the sequence is stationary and ergodic, 
from the time evolution of y^ m '(t) the probabilities 

P({We }) are computed, and one defines the block en- 
tropies of size e: 

H m (e) = - P { W c m) ) ^ (W M )) • (3) 

{W e (m) } 

Finally one introduces the e-entropy h{e): 

h(e) = lim h m (e) (4) 

rn — >oo 

where h m (e) = H m+ i(e) — H m (e) represents the e-block 
entropy growth at word length m. In a rigorous ap- 
proach, all partitions into elements of size smaller than 
e should be taken into account, and then h(e) is defined 
as the infimum over all these partitions |20|. The KS 
entropy can be identified as the limit e — > 0: 

h K s = \imh{e). (5) 

c->0 

In a deterministic chaotic system, one has hxs < 00 > 
in a regular motion h^s = 0, while for a random process 
with continuous states /i^s = oo. For some stochastic 
processes, it is possible to give an explicit expression for 
h(e) For instance, for a stationary Gaussian process 
with spectrum S(uj) ~ uj~( 1+2a ) with < a < 1 one has 

h(e) ~ e" 1 /" (6) 

while for i.i.d. variables whose pdf is continuous in a 
bounded domain (e.g. independently distributed vari- 
ables in [0, 1]) one has: 

ft( e )~lnQ) (7) 

Of course, letting aside the problem of the periodicity 
induced by the discrete nature of the states, a PRNG 
is good when its hxs is very large, such that the un- 
certainty on the "next" outcome is larger and the deter- 
ministic constraints appear on scales smaller than an e c 
defined by: {e c ) d ~ e~ hKS . 

On the other hand, in data analysis, the space where 
the state vector x lives is unknown and typically in exper- 
iments only a scalar variable u(t) is measured. Therefore, 
in order to reconstruct the original phase space, one uses 
the vector 

y< ra) (i) = (it(t), u(t + 1), . . . , u(t + to — 1)) G R™ 1 ; (8) 



that is another way to coarse grain the phase space. In 
this case, i.e. looking only at one variable, the maxi- 
mum scale where the Kolmogorov-Sinai entropy may be 
revealed is given by Cgi ~ e~ hKS , that is much smaller 
than e c , for large d 9j. Moreover, the single- variable 
words length that is necessary to consider, in order to 
detect hpcSi must be greater than d. This effect is harm- 
ful from the perspective of data analysis, but is welcome 
here. 

We also note that in any series of finite length T, it 
is not possible to have a good statistics of TO-words at 
resolution e if T < e mh ( e \ Therefore, for almost all the 
practical aims, i.e. for finite e and finite size of the se- 
quence, a chaotic PRNG with very high hxs has entropic 
properties indistinguishable from those of a perfect RNG. 

A. High entropy PRNG 

Several PRNGs are indeed discrete state versions of 
high entropic dynamical systems. A popular example is 
the multiplicative congruential method: 

z(t + 1) = az(t) mod M (9) 

where z(t), a and M are integers, with M > a > 1. 
In the following the term "map" will denote a dynam- 
ical system with continuous state space, and the term 
"automaton" a system with discrete state space (we al- 
ways assume a discrete time) . To avoid confusion, we will 
use the symbols z(t),w(t),z(t), w(£) only for discrete dy- 
namical variables (also vectorial) and x(t),y(t), x(t),y(t) 
for real dynamical variables. Eq.@ corresponds to the 
chaotic map: 

x(t + 1) = ax(t) mod 1 (10) 

where x(t) = z(t)/M. It is easy to see that the system 
has a uniform invariant p.d.f. in [0, 1] and has — 
In a. Therefore, only looking at e < e c ~ 1/a, one can 
capture the deterministic nature of the PRNG . 

It is worthwhile to stress that the chaotic features of 
the automaton are apparent only if one observes the sys- 
tem (0 after a coarse-graining procedure, namely with 
e ^> 1/M. Below this level of observation, the system 
keeps a loose trace of the chaotic features of its continu- 
ous precursor and, at the maximal resolution, at the first 
time step the block entropy already assumes its maxi- 
mum value H m (l/M) « InT for all to > 1 indepen- 
dently on the value of hxs- This happens because we 
are observing the complete state of the (one-dimensional) 
automaton and suggests again that a suitable attitude is 
to extract partial information from high-dimensional sys- 
tems. 

In the next subsection we show an alternative way, 
based on the high-dimension effect, to produce random 
number (up to a given word length) with an automaton, 
even at the finest resolution achievable. 
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B. High dimension PRNG 

It is known that non-chaotic high dimensional systems 
may display a long irregular regime as a transient effect 
[l4| . In this subsection, we show that, with a proper 
use of a transient irregular behavior, also systems with a 
moderate has may successfully generate pseudo-random 
sequences: in these cases, one observes a transient in the 
block e-entropies H n (e), characterized by a maximal (or 
almost maximal) value of the slope H n (t)/n, and then a 
crossover to a regime with the slope of the true hxs of 
the system. 

The most used class of PRNGs using this property 
are the so-called lagged Fibonacci generators 15], which 
correspond to the following map: 



x(t) = ax(t — n) + bx(t — T2) 



mod 1 



where a and b are 0(1) and t± <t 2 . 

Notice that Eq.JTTJ can be written in the form 

y(f) - Fy(t - 1) 

where F is a t 2 x t 2 matrix of the form 



(11) 



(12) 





a 
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F = 





1 







V 



(13) 



1 J 



showing explicitly that the phase space of l|ll|) has dimen- 
sion t 2 . It is easy to prove that this system is chaotic 
for each value of a, b £ N, with a, b > 0. The KS- 
entropy does not depend on t\ and T2 and is of the order 
~ ln(a6); this means that to obtain high values of hxs 
we are forced to use large values of o, b; nevertheless, the 
lagged Fibonacci generators are used with a = b = 1 . For 
these values of the parameters e~ ,lKS w 0.618 and e c i is 
not small. This implies that the determinism of the sys- 
tem should be detectable also with a large graining. De- 
spite these considerations, these generators work rather 
well: the reason is that the m-words, built up by a single 
variable of the T2-dimensional system (|12|) . have the 
maximal allowed block-entropy, H m {e) = mln(l/e), for 
m < t 2 , so that: 



H m (e) 



mln(-l) for 
r 2 In (i) + h KS (m - r 2 ) for 



m < T2 
m > T2 



(14) 



Ea. (|14f> has the following interpretation: though the 
"true" fiKs is small, it can be computed only for very 
large value of m. Indeed, by observing the one- variable 
m-words, which corresponds to an embedding procedure, 
before capturing the dynamical entropy one has to real- 
ize that the system has dimension r 2 , and this happens 
only for words longer than r-x- Fig. ^ shows H m (e) for 
T\ = 2, T2 = 5 and different values of e. 
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FIG. 1: The e block-entropy for the Fibonacci map with 
Ti = 2, rz = 5 , a = b = 1 and different values of e. The 
change of the slope from ln(l/e) to h,KS is clearly visible for 
m — T2 = 5. 



The importance of the transient behavior of H m has 
been underlined by Grassberger |22j who proposed an- 
other quantity beyond the KS entropy: the "effective 
measure of complexity" , namely 



C = 



=1 



m(h ri 



h m )- 



(15) 



From the above definiton, it follows that for large m, the 
block entropies grow as: 



C + mh 



KS- 



(16) 



For trivial processes, e.g. for Bernoulli schemes or 
Markov chain of order 1, C = and has > 0, while 
in a periodic sequence hxs — and C ~ bi(T). In the 
case of Fibonacci map, for small e, 



C = t 2 



In 



hKS 



t 2 In 



(17) 



For large r 2 (usually values O(10 2 ) are used) C is so huge 
that only an extremely long sequence of the order exp(r2) 
(likely outside the capabilities of modern computers) may 
reveal that that the "true" KS entropy is small. 

Let us now discuss the behavior of the discrete Fi- 
bonacci generator: 



z(t) =az(t-T 1 )+bz(t-T 2 ) mod M (18) 

where z(t) £ [0, M — 1] and M ^> t 2 - The parameters t\, 
T2 and M are chosen in order to have a period as long 
as possible. Number-theoretical arguments 0] allow to 
choose these parameters such that the period of the orbit 
is maximum T = M T2 — 1 . 

When the period is maximum, for e > l/M one has: 



rn 



H m (e) 



t 2 In (i) + h KS {m - t 2 ) for 
r 2 ln(M) for 



for m < 



T 2 

t 2 < m < m* 



m > m* 



(19) 
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where 



't-ks 



In [ - 



InM + h 



KS 



(20) 



When e = 1/M we have m* = ra, the second regime 
in Ea. H19|) disappears and the block entropy behavior is 
independent of hxs- Still, as for the continuous case, if 
T2 is large one observes only the pseudo-chaotic transient 



is only quadratic in the deviation (this is a simple con- 
sequence of the fact that the value of the entropy in Eq. 
(|22|l is the maximum achievable). Thus, it is really diffi- 
cult to numerically observe unbalances in the frequency 
of some specific words by studying only the block en- 
tropies. 



A. The spectral test 



fr m (e)«m]nl-J. (21) 

Summarizing, systems with high values of hxs or with 
high dimension, produce sequences having entropic prop- 
erties rather close to those of a perfect RNG, for two 
different reasons. In the first case, the large Iiks allows 
to use a small e (but large enough to achieve a proper 
coarse-graining): in such a way, the e-entropy coincides 
with that of a perfect RNG. In the second case, the high 
dimensionality of the system prevents the entropic anal- 
ysis to reveal the asymptotic value of h^s before the 
end of a long transient behavior mimicking a complete 
random system. 

In conclusion, a deterministic PRNG has good entropic 
properties for long (but finite) sequences if Iiks or C are 
large. In Section Hvl we will propose a multi-dimensional 
cat map as a PRNG having both these properties. 

III. TESTS FOR PRNGS 

Several techniques have been developed in order to test 
"how random" is a given sequence of numbers. These 
algorithms are available in easy-to-use software pack- 
ages collecting dozens of different tests, like, for ex- 
ample, the DieHard [23| and the NIST batteries. 
Many of them compute the frequency of some words 
f(z(j), z(j + 1), . . . , z(j + n)) made up of n consecutive 
outputs of the generator and compare it to the theoret- 
ical probability in the random case. Examples are: the 
frequency and block- frequency tests (computing f(z(j)), 
i.e. m = 1), Poker test looking for words with to = 5 cor- 
responding to the Poker hands (e.g. Full house: 00011, 
four a kind: 00001), template tests checking the occur- 
rences of some (~ 10 2 ) words with length m = 10 — 12. 

It is worth stressing that all these benchmarks are au- 
tomatically passed if the block e-entropy, is maximal for 
words of length to, namely 

H m (e) =mln(M). (22) 

with e = 1/M where M is the number of symbols pro- 
duced by the generator. When H m is maximal, it is im- 
possible to distinguish the output from a truly random 
sequence by looking only at to consecutive symbols. The 
reason for introducing many tests, instead of looking only 
at the entropy, is that a little departure from a constant 
word frequency implies a correction in the entropy that 



The entropic analysis is a very powerful tool from a 
theoretical point of view but it presents a major limi- 
tation: in a phase space with a finite number of states 
the block entropy cannot grow more than In T, where T 
is the period of the orbit. This fact essentially fixes an 
upper bound on the number of possible equidistributed 
words and on their length. Even with the longest pe- 
riods available in current used algorithm, the bound on 
the length is of order 10 2 — 10 3 . On the other hand, 
computer simulations often necessitate of a large amount 
of random numbers, and long-term relationships among 
these numbers can be sources of hard-to-discover biases 
[26l I2H |28| . Therefore, less severe tests than the en- 
tropic one are needed. One can ask that the correlations 
among different outputs (or, more generally, among dif- 
ferent functions of the outputs), vanish even when the 
outputs are at a time distance greater than the scale 
where Ea. H22|) ensures the equidistribution of the words. 
On this timescale, numbers should appear random, as 
far as one is interested on statistical observable, even if, 
looking at the whole sequence, later numbers are com- 
pletely determined by previous ones. Indeed, according 
to the entropic analysis, the knowledge of approximately 
InT consecutive outputs permits to determine exactly 
the discrete starting condition of the system and con- 
sequently to predict the whole sequence, removing any 
randomness from it. 

The main tool to analyze the property of correlation 
functions is the spectral test. We start defining the fre- 
quency f{z(ti), . . . , z(t n )) of the word (z(ti), . . . , z{t n )) 
as we made for the Kolmogorov entropy, where now tj's 
are generic times and z(ij)'s are not in general consecu- 
tive outputs. The spectral test is the multi-dimensional 
Fourier transform of / (z(ti), . . . , z(t n )) 

/0i,S2, ,s n ) = ^exp {^^s 3 z{t 3 )^ ^ , (23) 

where Sj £ [0, M — 1], M is the number of the discrete 
states and (...) denote the average over the trajectory. 
For true random numbers: 

/(si, s 2 , , s n ) = S slfi S S2! o ■ ■ ■ <5s„,o (24) 

for any choice of n and of the time lags Vs. Values of 

the function /(si,s 2 , ,s n ) significantly different from 

denote wave vectors of probability density fluctuations 
in the lattice z(ti), z{t2) ■ ■ ■ z(t n ). These fluctuations can 
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be safely neglected only when their characteristic length 
scale (which we assume to be the inverse of the modulus 
of the wave vector) is much smaller than the maximum 
precision one is interested in. Many generators (i.e. the 
linear congruential class of gene rators) produce numbers 
that "fall mainly in planes" [171 flSl , and the presence 
of these planes is detected by the spectral test. 

The importance of the spectral test is related to the 
fact that analytical or semi-analytical methods 0, |2jJ 
allow for a fast calculation for simple systems. Further- 
more, since any L 2 function can be written as a Fourier 
series, condition (|24|l implies the vanishing of any corre- 
lation of up to n functions of time-delayed variables: 



(gi(z(t 1 ))g 2 (z(t 2 ))...g n (z(t n ))) = 
(g 1 (z(t 1 )))(g 2 (z(t 2 ))) . . . (g n (z(t n ))} 



(25) 



for every <?.; 6 L 2 



IV. THE CAT-MAP AS A RANDOM NUMBER 
GENERATOR 



Ea. (l28|) is symplectic, indeed one can write Ea. (|28() and 
Ea. (|27|) as a canonical transformation 



x= aS(x',y) , = dS(x!,y) 
dy ' 9x' 

where the generating function is given by 



(29) 



N N 

<?(x',y)=5>^-~ ^ 



2 ^ 



(y J A jk y k +x' j B jk x' k ). (30) 



It can be shown that when Tr(M) > 2N map (O is 
an Anosov system with uniform invariant measure. The 
output of our generator will be the first component of the 
vector x. The condition N > 1 raises the Kolmogorov 
entropy, cancels the correlation and increases the length 
of the periodic orbits (in the discretized case). We will 
describe in detail these three aspects in the following. 

First of all, according to Pcsin identity, the Kol- 
mogorov entropy of this system is equal to the sum of 
the positive Lyapunov exponents: 



Recently some authors [3(J proposed the use of the 
Arnol'd cat map as a PRNG. We will briefly recall 
the properties of this map and then propose a multi- 
dimensional version with N coupled maps, showing that 
this generalization gives rise to very good statistical prop- 
erties. In particular, it has both the properties analyzed 
in Section ^ giving maximal e-entropy, namely it pos- 
sesses a high value of Iiks and it is a high-dimensional 
system. We will see that this system has also very good 
properties from the point of view of correlation functions. 

The 2-dimensional Arnol'd's cat map [3l| is a sym- 
plectic automorphism on a torus satisfying the property 
of Anosov systems 32], namely it is everywhere hyper- 
bolic and has a positive Kolmogorov entropy. The map 
reads 



IKS 



1 

b 1 



ab 



mod 1, 



(26) 



where a, b 6 N. The standard example given by Arnol'd 
is obtained with a = b = 1. 

The multi-dimensional generalization can be written 
in the following way: 



mod 1, 



with 



A 



(27) 



(28) 



where M is a 2N x 2N matrix, x, y e l", I is the N x N 
identity matrix and A, B are symmetric N x N matri- 
ces with natural entries in order to obtain a continuous 
mapping. It easy to see that the evolution law given by 



Ai>0 



A; 



(31) 



A d-dimensional hyperbolic symplectic system possesses 
exactly d/2 positive Lyapunov exponents; if they are of 
the same order of magnitude, the Kolmogorov entropy 
grows proportional to the number of dimensions. Notice 
that the Kolmogorov entropy may also be raised by sim- 
ply taking the matrix A, B with very large entries. This 
method, however, produces only an increase in the en- 
tropy which is logarithmic in the size of the entries. Of 
course for the system 1)271- 125)1 the A, are easily obtained 
from the eigenvalues on of M: Ai = In | cci | . 

For the 2D cat map we compute an approximate value 
of the e-entropy, obtained as H^(e) — H^(e) varying e 
and the parameters a, b. In order to highlight the prac- 
tical limitations of PRNG we study the discrete version 
of Eq. 11261) with M = 2 20 possible values of x and y (see 
the next section for details). 

As one can observe in Fig|3 both the standard prob- 
lems of PRNGs appear. Indeed the figure shows that 
decreasing the value of e we observe a "plateau" around 
the value of hpes- At lower values of e, there is an abrupt 
decrease due to the periodic nature of the map. Never- 
theless it seems that if we use the map as a generator of 
a number of symbols ~ 1/e with e > e c we are, with a 
good approximation, near the value corresponding to a 
theoretical RNG, given by h(e) = — ln(e). Let us note 
that, when hxs is large enough (curves for a = 5, b = 7, 
Iiks ~ 3.61 and a = 11, b = 17, Hks — 5.24), because 
of the limited number of allowed states, one does not 
observe the plateau h(e) ~ h^s- 

Let us study the properties of the time correlation 
of the outputs. The following result holds: Let ei be 
the 2N- dimensional vector (1,0,0...). If the vectors 
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FIG. 2: H 4 (e) - H 3 (e) for the 2D cat map of Eq.lO with 
M = 2 20 , different values of a, 6 as a function of e. The 
horizontal lines indicate the has values 

(M T ) tl e 1 , (M T ) t2 ei, . . . , (M T )* 2N ei are linearly indepen- 
dent, then one has: 

exp ^27ri^ sjxxitj^j ^ =S Sli0 S S2 fl---S S2Ni0 . (32) 

Furthermore, the independence of the vectors is ensured 
for any choice of the time delays ti 's if the matrix M 
has real positive and non- degenerate eigenvalues and the 
vector ei has non-zero component on all the eigenvectors. 
For the proof see Appendix A. 

The practical meaning of this result is the following: 
we observe only the variable Xx, keeping the remaining 
2N — 1 variables hidden, and we study its correlation 
functions. In this way correlation functions involving up 
to 2N different times vanish, i.e. Ea. (|25[) holds for n < 
2N, because the contributions due to different values of 
the hidden variables cancel out in the averaging. 

The result of Ea. (|32[l can be taken as one of the 
strongest characterization of a finite random sequence: 
as we said, word equidistribution can hold only up to 
a value of n ~ InT where T is the length of the se- 
quence. Indeed, some authors |25l ] define as random se- 
quence of length T one containing all the possible words 
up to length n. On the other hand, Eq.JH^J) is a gen- 
eralization of that condition: for consecutive time delays 
tj = j the two properties are equivalent, while for generic 
values of the tj's it ensures long-range independence of 
the outputs, without asking an exponential number of 
equiprobable words in the sequence. 

The validity of the property of (|3*2|l in the discrete case 
will be subject of careful analysis in the following section. 
Here, we just point out that, even in the continuous case, 
this property is not shared by some of the dynamical sys- 
tems used for generating random numbers. For example, 
let us recall the Fibonacci map 

x(t) — ax(t — Ti) + bx(t — t-i) mod 1 (33) 



with a, b 6 N with a, b > and t 2 > T\. It is straightfor- 
ward to show that the correlation function 

(exp{2ni(six(t) + s 2 x(t - n) + s 3 x{t - r 2 ))}) (34) 

is not equal to zero for the vector s = (si,S2,S3) = 
(1, — a,— b). Thus, even if the dimension of the phase 
space Ti may be very high, it is sufficient the three-points 
correlation function to unveil the deterministic nature of 
the system. 

Therefore, when the dimension of the cat map, 2N, 
is equal to the dimension of the Fibonacci generator t 2 , 
both the systems guarantee that words of length 2N are 
equidistributed. However, the main advantage of the 
multi-dimensional cat map is that also the words made up 
of 2N non-consecutive symbols are equidistributed. This 
property does not hold for Fibonacci generators and in 
some case this can lead to serious problems. A famous 
example is the "Ferrenberg affair" [2||: persistence in 
binary Fibonacci generators gives misleading results in 
Montecarlo simulation. This problem is well analyzed in 
the framework of information theory in |27| . 

In the next section we numerically study the discrete 
version of the multidimensional cat map and we check 
whether the good statistical properties of the system sur- 
vive in this case. 



V. NUMERICAL ANALYSIS AND TEST OF 
THE MULTIDIMENSIONAL CAT AUTOMATON 

A digital computer cannot handle real numbers. What 
a computer really calculates is a finite-digit dynamics 
that can be represented as a dynamics on integers. We 
will consider in the following the multidimensional cat 
automaton, namely 

(w') = (i i+ A ba)(w) modM ' ( 35 ) 

where, as usual, A, B have natural entries, Zi,wi 6 
[0,1,...,M-1]. 

The first problem in passing from the continuous to 
the discrete case is that the system has a finite number 
of state M 2N and consequently it must be periodic and 
it is no more truly chaotic. The optimal condition is that 
there is only one orbit covering all the states but the ori- 
gin z = w = (which is a fixed point) thus obtaining a 
period T = M 2N — 1. A peculiar feature of cat maps 
is that periodic orbits of the continuum system have ra- 
tional coordinates |3^| . Consequently, the orbits of the 
discretized version of the map are completely equivalent 
to periodic orbits of the continuous system with coor- 
dinates Zi/M,Wi/M . As a corollary, being the map in- 
vertible, periodic orbits do not have any transient: every 
state is recurrent. 

Unfortunately the cat map has typically many orbits, 
and the great majority of them are of the same length. 
A theoretical analysis of these orbits have been made 
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[33j in the 2-dimensional case. For generic dynamical 
system, arguments based on random maps suggest that 
the average length of the orbits T should be roughly the 
square root of the total number of states [j| : in our case 
T « M . This scaling have been numerically observed 
in typical chaotic dynamical systems j^il]. in Fig- 01 
we show the length, T, of the orbits as a function of 
M and N. Despite the large fluctuations, one retrieve 
the expected qualitative scaling and, more importantly, 
the periods seem to be independent on the initial condi- 
tion; this suggests that several symmetry operations exist 
mapping one orbit into another one of the same length. 
Since the choice of M is critical in determining the length 
of the orbits, we restrict ourselves to prime number values 
of M , in order to avoid the presence of trivial invariant 
sub-lattices generated by the divisors of M. Notice also 
that data are plotted as a function of M N : the lengths 
show the correct exponential growth with N , at least in 
a statistical sense. We stress that in Fig. |3we show the 
result for M not very large (< 10 3 ), from the observed 
behavior one can say that for M ~ 10 9 the period should 
be extremely large. 
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FIG. 3: The period of the orbits for the multidimensional cat 
map, with different values of N, as a function of M (only the 
prime values of M are considered). The straight lines refer to 
the probabilistic argument T ~ M . 

Unfortunately we have no theoretical control over the 
period, and wild fluctuations are present when M varies; 
therefore it is better to choose a value of M, N, A, B 
and directly check the value of T or a lower bound. 
In the following we will consider the choice N = 3, 
M = 1,001,400,791 and 



(36) 



With these parameters, the hypothesis leading to Eq. H32|) 
hold, furthermore we numerically obtained T > 7 • 10 12 , 
that is a satisfying lower bound for typical simulations. 

The very encouraging result we obtained for the cor- 
relation functions in the continuous case is the main rea- 
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son for the use of the multidimensional cat automaton 
as a PRNG, but we have to check whether the property 
proven in the previous section holds also in the discrete 
case. The choice of <|36[) satisfies the hypothesis of the 
theorem, namely the eigenvalues are all strictly positive, 
non-degenerate and the vector (1,0,0,0,0,0) has non- 
zero component along all eigenvectors. 

From a general point of view there are two main tech- 
nical caveats when passing from continuous to discrete 
systems concerning the theoretical spectral test. In our 
case, since the orbit do not cover all the space, it is a pri- 
ori impossible to average on the uniform distribution as 
discussed in the Appendix for the continuous state case. 
Even if we suppose that the orbit is sufficiently homoge- 
neous that Ea. fiUI) (see the Appendix) is approximately 
valid, condition (|42J) becomes in the discrete case 



2N 

£■ 



l lk 



mod M Vfc = 1, . 



,27V 



(37) 



since it is sufficient to consider vectors (si, S2, ■ ■ ■ , S2n) 
in the first Brillouin zone. 

These two reasons prevent us to perform the spectral 
test using the nice theoretical arguments used for other 
kind of generators |19L l35j and force us to use numerical 
simulation. This constitutes a hard computational task 
and we perform the test for low values of M and studying 
products of the form 

/(si,s 2 ) = <^exp (~^ s ^( t ^)] ex P (~^ S2Zl ^ 

'(38) 

We use M = 1031 and N = 3 obtaining a period 
T = 274,243,921 and letting Sl ,s 2 G [0, M - 1]. Us- 
ing an FFT numerical algorithm we check up to time 
delays ti—t\ < 250 that for all values, but si = s 2 = 0, 
|/(si,S2)| < 10~ 5 . With a lower number of states, 
M = 127, T = 1,016,188, we compute also the three- 
point spectral test, obtaining always values compatible 
with the inverse of the square root of the period T. 
This suggests that the periodic orbits look like a fi- 
nite statistical sample of the continuum equilibrium dis- 
tribution, as far as one studies only few-point correla- 
tion functions. An important remark is that a solution 
(si, s%, . . . , Sjjv) °f t nc Diophantine Ea. (|37J) implies that 
/(«!, s|, . . . , sj-jiv) = 1 independently on the stationary 
distribution and, consequently, on the period T. This 
means that the low values observed in the numerical spec- 
tral test exclude the possibility of a solution in Eq. I|37l) . 

In order to look for any other possible bias, we also 
apply the NIST battery to test our multidimensional cat 
automaton, with the parameters of Ea. (|36|l . for generat- 
ing 10 3 binary string of 0's and l's of length 10 6 ; all tests 
performed with the recommended parameters have been 
passed. 
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VI. CONCLUSIONS 



APPENDIX A 



In this paper we show how, using properties of high 
dimensional deterministic chaotic systems, it is possible 
to generate a good approximation of a random sequence, 
in spite of unavoidable constraints of deterministic algo- 
rithms running on digital computers. 

Summarising, we have two possible mechanisms to ob- 
tain good PRNGs using deterministic systems: very high 
KS entropy, and "transient chaos" with a large finite-time 
e entropy, due to the high dimensionality of the system. 
We propose the multi-dimensional cat map as a PRNG 
having both these properties. Another important exam- 
ple of system with both the properties is the one proposed 
by Knuth 36] : one iterates the Fibonacci generator 
with M = 2 31 — 1, n = 37 and r 2 = 100, with this choice 
the period is extremely large, then the output sequence 
is obtained taking the variable in Ea. (|llf) every T steps 
(T = 1009 or 2009). In such a way, for the words of size 
up to T2 (i.e. extremely huge), the e entropy is practically 
~ ln(l/e), i.e. like a perfect RNG. Moreover, even if the 
simple Fibonacci generator fails the three-points spectral 
test, it is harder to find non-null vector in the spectral 
test of Knuth's generator, because of the fact that T, 
Ti and T2 are relatively prime numbers. Nevertheless, it 
does not seem that a general result like that of Ea. H32|) 
may be easily extended to this PRNG. 

We suggest that the multi-dimensional cat map is suit- 
able for generating random number sequence. The main 
advantage if compared with other generators, is the fac- 
torization of all n-times, with n < 2N, correlation func- 
tion due to the high dimensionality of the system and 
the presence of hidden variables. This result is rigorously 
true (also in the case n = 2N) in the continuous system; 
numerical checks show that this property survives in the 
discrete case. Moreover, this map has a large value of the 
KS entropy giving good entropic properties at non-zero, 
but small, e. 

A disadvantage of this method is that we can not 
predict analytically the period given the parameters or, 
equivalently, write a condition on the parameters in or- 
der to obtain the maximum period. However, probabilis- 
tic arguments |16| . confirmed by numerical check, show 
that the period increases exponentially with N, there- 
fore with a proper choice of the parameters we achieve 
extremely large periods. An analytical criterium to pre- 
dict the length of the period could pave the way to the 
application of multi dimensional cat maps as high quality 
PRNGs. 



Consider the system of Eqs. (|27() and (|28() . In this 
Appendix we give the proof of the following proposition: 

Let ei be the 2N- dimensional vector (1,0,0...). // 
the vectors (M T ) tl e 1 , (M T )* 2 ei, . . . , (M T )* 2Jv ei are lin- 
early independent, then one has: 



2N 

exp | 2iri sjxi (tj) 



= 6 



s 1 M°s 2 ,0 ■ ■ 



■$S 2 



(39) 



Furthermore, the independence of the vectors is ensured 
for any choice of the time delays ti if the matrix M has 
real positive and non-degenerate eigenvalues and the vec- 
tor ei has non-zero component on all the eigenvectors. 

Proof. Since the system under study is ergodic and its 
invariant measure is uniform, we can write the average 
in Ea. l|32)l as: 



27V 



dx\ . . . dx2N exp 2iri s j x i(tj) 



(40) 



Here, with an abuse of notation, we call Xj the compo- 
nents of both the x and the y vectors, i.e. Xj+N = Vj, 
, N. Let us also call m^j. the elements of the 
: . We can rewrite the previous expression in 



1. 



matrix 



the following way: 



dx\ . . . dx2N exp 




lk 



(41) 



Notice that we do not take care about the modulus since 
the Sj are integers and the complex exponential is peri- 
odic. When integrating over the x^s, the result is zero 
for every value of the s^'s, excluding the values that are 
solution of the linear system: 



2N 

s^mf^ = 



VTc = 1 . . . 2N 



(42) 



that yield 1 as a result of Ea. (|41|l . Since sj = Vj 
is a trivial solution for the linear system, it is sufficient 
to show that this solution is unique to demonstrate Eq. 
(|32|l . In particular, by Cramer's rule, it is sufficient to 
show that det G^O, where G is the matrix of coefficients 
of the linear system (|42|l . namely 



9kj 



{h) 
= m ik • 



(43) 
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Notice that the column of the matrix G are constituted by 
the component of the vectors (M T ) tj ei; this means that 
the condition of having det G ^ is equivalent to require 
that the vectors (M T )' 3 ei are linearly independent. This 
completes the first part of the proof. 

Now, we show how Eq. (|42|) has always a unique solu- 
tion when the matrix M has positive and non-degenerate 
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eigenvalues A&, and the vector ei has non-zero compo- 
nent along all the eigenvectors, In this case, we rewrite 
the matrix G in the eigenvector basis, obtaining 



/ A 



det G = C\C2 ■ ■ ■ C2N det 



4 1 



A? 



2 A 



A* 2JV \ 
A* 2W 



\*2JV , 

A 2JV / 



(44) 



where Ck ^ for hypothesis is the component of ei along 
the k-th eigenvector. Notice that the c^'s are real: since 
the eigenvalues are real by hypothesis, also the eigenvec- 
tors have real components in the natural basis of M. 2N . 
The proof is a reduction ad absurdum: let us suppose that 
det G = 0, this implies that there exist a linear combina- 
tion of the columns satisfying 



E 



W = 



Vk = 1 . . . 2N. 



(45) 



Calling P(z) the polynomial 



(46) 



Ea. (|45|) implies that the polynomial P(z) has 2N distinct 
real positive roots being, by hypothesis, all eigenvalues 
positive and non-degenerate. Then, by Descartes sign 
rule, it must have at least 2N sign changes in the coef- 
ficients, but this is impossible, since P{z) has just 2N 
terms different from 0. Thus, det G is necessarily differ- 
ent from zero for any possible choice of the time delays; 
this completes the proof. 
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